#!/bin/csh
# 1d_timeaverage.csh
# After producing a file with the time-series of the vertical profiles, it calculates 
# the time-averaged profiles over the time interval time1-time2, The resulting file
# is timeave_profiles.dat is moved back to the parent simulation directory.
#
# Run it as ./1d_timeaverage.csh time1 time2
# where time1 and time 2 are the initial and final non-dimensional times 
# over which the time-averages must be computed
##################################################################################

#-------------------------------------
# Compute the profiles
#-------------------------------------
source simulation.dir

set tzfile = time_profiles.dat

set time1 = ${1}
set time2 = ${2}

# Check whether 'time_profiles.dat' already exists and, if not, create it using make_timeprofiles.csh
if (! -e {$dirout}$tzfile) then
  echo 'Run make_timeprofiles.csh'
  ./make_timeprofiles.csh 
endif

# Determine from the first profiles.dat file the number nz of vertical points
set file = (`ls -ltr ${dirout}profiles_* | head -1| awk '{print $NF}'`)
set nz = `awk 'END { print NR-1 }' $file`

# Compile the fortran code used to actually compute the time averages...
echo 'Compile fortran code...'
gfortran -o timeave 1d_timeaverage.f90
#... and run it to generate the file timeave_profiles.dat
echo '...and run it'
./timeave $dirout $time1 $time2 $nz

# Move timeave_profiles.dat into the parent simulation directory
mv timeave_profiles.dat $dirout 


#------------------------------------------
# Plot the time-averaged profiles with GMT
#------------------------------------------
set dirplots0 = {$dirout}Plots/
set dirplots = {$dirout}Plots/Timeaveraged_profiles/ 

set makejpg = 1
set jpgdensity = 150
set deletepdf = 1

#---------------------------
# Select profiles to plot
#---------------------------
set temperature = 0
set velocity = 0
set urms = 1
set ux = 1
set uz = 1
set viscosity = 0
set composition = 1

if (! -d $dirplots0) then
  mkdir $dirplots0
endif

if (! -d $dirplots) then
  mkdir $dirplots
endif

set profiles = ${dirout}timeave_profiles.dat
echo "Plot:" $profiles

#----------------------------
# Set plot size in cm
#---------------------------
set dimx = 10
set dimz = 10

#----------------------------
# Set plot origin (in cm)
#----------------------------
set X0 = 3
set Y0 = 5

#--------------
# Set ranges 
#--------------
set zmin = 0
set zmax = 1

#-----------------------------------
# Set temporary ps output file
#------------------------------------
set output = timeaveprofile.ps

#----------------------------------
# Set gmtdefaults
# ---------------------------------
source setgmtdefaults

####################
# Temperature
####################
if ($temperature == 1) then
  echo "Time-averaged temperature profile" 

  set title = ""
  set xlabel = "Temperature"
  set zlabel = "Depth"

  set pdfoutput = temperature_timeaveprofile.pdf
  set jpgoutput = temperature_timeaveprofile.jpg

  set min = `awk '{if(min==""){min=$2}; if($2< min) {min=$2}} END {printf "%12.5f", min}' $profiles`
  set max = `awk '{if(min==""){min=max=$2}; if($2>max) {max=$2}} END {printf "%12.5f", max}' $profiles`
  set delta = `echo "scale=2; ($max - $min)/5" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba${delta}:${xlabel}:/a${deltaz}:${zlabel}:WSen:.${title}: -K >! $output
  awk '{print $3, 1-$1}' $profiles  | \
  psxy -R -JX -W7/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

####################
# Velocity
####################
if ($velocity == 1) then
  echo "Velocity profile" 

  set title = ""
  set xlabel = "Velocity"
  set zlabel = "Depth"

  set pdfoutput = velocity__timeaveprofile.pdf
  set jpgoutput = velocity__timeaveprofile.jpg

  set min = 0
  set max = `awk 'NR > 1 {if(min==""){min=max=$2}; if($2>max) {max=$2}} END {printf "%12.5f", max}' $profiles`
  set deltax = `echo "scale=0; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba{$deltax}:${xlabel}:/a${deltaz}:${zlabel}:WSen:.${title}: -K >! $output
  # RMS 
  if ($urms == 1) then
    awk 'NR>1 {print $2, 1-$1}' $profiles| \
    psxy -R -JX -W7/0/0/0 -O -K >> $output
  endif
  # Horizontal 
  if ($ux == 1) then
    awk 'NR>1 {print $7, 1-$1}' $profiles| \
    psxy -R -JX -W7,0/0/255 -O -K >> $output
  endif
  # Vertical 
  if ($uz == 1) then
    awk 'NR>1 {print $8, 1-$1}' $profiles| \
    psxy -R -JX -W7,255/0/0 -O -K >> $output
  endif

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif


####################
# Viscosity
####################
if ($viscosity == 1) then
  echo "Viscosity profile" 

  set title = ""
  set xlabel = "log@-10@-(Viscosity)"
  set zlabel = "Depth"

  set pdfoutput = viscosity_timeaveprofile.pdf
  set jpgoutput = viscosity_timeaveprofile.jpg

  set min = `awk 'NR > 1 {if(min==""){min=log($4)/log(10)}; if(log($4)/log(10)< min) {min=log($4)/log(10)}} END {printf "%12.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=log($4)/log(10)}; if(log($4)/log(10)>max) {max=log($4)/log(10)}} END {printf "%12.5f", max}' $profiles`
  if ($min == $max) then
    set min = -1
    set max = 1
  endif
  set deltax = `echo "scale=1; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba{$deltax}:{$xlabel}:/a${deltaz}:{$zlabel}:WSen:.{$title}: -K >! $output
  awk 'NR>1 {print log($4)/log(10), 1-$1}' $profiles | \
  psxy -R -JX -W7/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

######################
# Composition
######################
if ($composition == 1) then
  echo "Composition profile" 

  set title = ""
  set xlabel = "Composition"
  set zlabel = "Depth"

  set pdfoutput = composition_timeaveprofile.pdf
  set jpgoutput = composition_timeaveprofile.jpg

  set min = `awk 'NR > 1 {if(min==""){min=$6}; if($6< min) {min=$6}} END {printf "%12.5f", min}' $profiles`
  set max = `awk 'NR > 1 {if(min==""){min=max=$6}; if($6>max) {max=$6}} END {printf "%12.5f", max}' $profiles`
  if ($min == $max) then
    set min = -1
    set max = 1
  endif
  set deltax = `echo "scale=2; ($max - $min)/4" | bc -l`
  set deltaz = 0.2

  psbasemap -X$X0 -Y$Y0 -R$min/$max/$zmin/$zmax/ -JX$dimx/-$dimz -Ba{$deltax}:{$xlabel}:/a${deltaz}:{$zlabel}:WSen:.{$title}: -K >! $output
  awk 'NR>1 {print $6, 1-$1}' $profiles | \
  psxy -R -JX -W7/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif


#############################
# Moving and cleaning
#############################

if ($makejpg == 1) then
  mv -f *.jpg $dirplots
endif

if ($makejpg == 1 && $deletepdf == 1) then
  rm -f *.pdf
else 
  mv -f *.pdf $dirplots
endif

rm -f *.ps 
